When not to avoid inbreeding: a gene’s eye view perspective

Author

Thomas Keaney, Arvid Agren and Hanna Kokko

Load packages

Code
# for tidy style coding and plotting

library(tidyverse) 
library(vroom) # to read lots of csv files at once

# more table options

library(pander) # for tables
library(kableExtra) # for scrolling tables
library(data.table) # for efficient handling of large dataframes

# making ggplot more powerful

library(MetBrewer) # for colour palettes based upon artwork housed at the MET
library(MoMAColors) # for colour palettes based upon artwork housed at MoMA
library(wesanderson) # for colour palettes based on wes anderson movies
library(PNWColors) # for colour palettes 
library(tidybayes) # for plotting distributions
library(stickylabeller) # labelling facets with strings in ggplot
library(geomtextpath) # for curved plot annotations
library(ggtext) # for markdown syntax in plot labels
library(patchwork) # for patching plots together
library(ggnewscale) # to reset scales in plots, allowing multiple fill arguments in ggplot

# for computation speed checks

library(profvis) # breakdown of complex functions
library(bench) # individual functions

The model

The seminal equation

To model the inclusive fitness gained from an inbred mating, three components that contribute to fitness are required:

  1. The number of offspring produced directly: \(n\)

  2. The reduction in fitness (number of offspring) due to inbreeding: \(-\delta n\)

  3. The indirect fitness gain (number of offspring) due to inbreeding: \(rn\), where \(r\) is the relatedness coefficient

Put together, the inclusive fitness from a single inbred mating is:

\[(1 + r)(1 - \delta)n\]

while fitness from a single outbred mating is simply \(n\).

When \((1 + r)(1 - \delta)n \gt n\) selection should favour a preference for inbreeding.

Solving the inequality for \(\delta\):

\[\delta \lt \frac{r}{1 + r}\] which for varying values of \(r\) looks like this:

Code
inbreeding_maximum_function <- function(r){
  r / (1 + r)}

parameters <- expand_grid(r = seq(from = 0, to = 1, by = 0.05),
                          delta = seq(from = 0, to = 1, by = 0.05))

r <- parameters %>% distinct(r)

inbreeding_equilibria <- 
  map_dfr(r, inbreeding_maximum_function) %>% 
  rename(depression_threshold = r) %>% 
  bind_cols(r)

inbreeding_equilibria %>% 
  ggplot(aes(x = r, y = depression_threshold)) +
  geom_line(linewidth = 0.8) + 
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = '_r_, the relatedness coefficient',
       y = ~delta~'(inbreeding depression)') +
  scale_x_continuous(expand = c(0, 0.009)) + 
  scale_y_continuous(expand = c(0, 0)) +
  theme_bw() +
  theme(text = element_text(size = 14),
        axis.title.x = element_markdown())

Code
# (prop fitness lost)\n that can be tolerated"

The parameter space above the curve shows where inbreeding avoidance should evolve, while the parameter space below the curve shows where inbreeding preference should evolve.

Accounting for sex differences in genetic architecture

As stated above, inclusive fitness in the absence of inbreeding depression is \((1 + r)n\). Here \(r\) represents the correlation between genotypes carried by interacting females and males, under the implicit assumption that loci appear equally in both sexes. However, given that there is sexual dimorphism in genetic architecture for many taxa, \(r\) does not sufficiently represent the correlation between genotypes for all loci.

To delineate differences in the effect of \(r\) for different regions of the genome, we multiply \(r\) with a new variable \(a\), the probability that a locus present in one sex is also present in the gametes produced by the other. Unlike \(r\) which is relative to the population mean relatedness, \(a\) is expressed as an absolute value ranging from 0 to 1.

The indirect component of fitness accrued by from an inbred mating becomes

\[ran\]

and inclusive fitness from an inbreeding event becomes

\[(1 + ra)(1 - \delta )n\]

Taking an allele found at a diploid autosomal locus as an example, all of the gametes produced by a relative possess this locus, where they could potentially carry alleles identical by descent. In this case \(a = 1\) and the indirect component of inclusive fitness is dictated solely by \(r\). The results for this autosomal scenario are presented in Parker (1979), Kokko and Ots (2006) and others who have explored this topic. In contrast, an inbreeding preference allele present at a locus on a Y or W chromosome has no opportunity to propagate any alleles identical by descent through inbreeding, as these chromosomes are not carried by the gametes of the opposite sex mating partner. In this case \(a = 0\). However, as inbreeding depression is a result of homozygosity for deleterious recessive alleles throughout the genome, the costs of inbreeding depression are born by all alleles carried by the individual. Conflict over the expression of inbreeding preference between alleles present on autosomes and those present on hemizygous sex chromosomes is immediately clear.

X- or Z-linked loci present an interesting intermediate case, with sex-specific values for \(a\). When the inbreeding locus is carried by the sex with homozygous sex chromosomes, \(a\) is half that of autosomal loci, whereas it does not depart from the autosomal case when the locus is found within the hemizygous sex. Using loci on the X as an example, those present in a XX female are only found in ~50% of a interacting males gametes, as the remaining 50% carry Y chromosomes (assuming an even primary sex ratio). When an X-linked locus is found in a male, an interacting female’s gametes all carry X chromosomes and \(a = 1\).

The X/Z situation is made additionally complex because there is an element of frequency dependence to the kin selected benefits. When an inbreeding allele on an autosome is rare, then the chance of a relative carrying two copies is low, whereas when the allele is common, this chance is much higher. Rarity therefore leads to similar fitness outcomes for autosomal and X/Z linked alleles (when present in the hemizygous sex), while commonality of the allele likely roughly equates to the conflicting situation outlined in the above paragraph. However, frequency dependence might not as relevant as I initially expected, because the inbreeding allele becomes very common quickly within families, even whilst rare across the population (first proposed in Fisher, 1930).

Table 1. Values of the parameter \(a\) for different regions of the genome. \(a\) is the one-way probability that a locus carried by one individual is found within the gametes of an opposite sex individual. Note that cytoplasmic chromosomes are assumed to have exclusive maternal inheritance.

Code
x <- 
  c(1, # autosomes, X chromosome males or Z chromosome females, haplodiploid both sexes when producing females 
    0, # Y or W chromosome
    0.5 # X chromosome females or Z chromosome males
  )

tibble(`Prob. that opposite sex gametes carry focal locus` = c(1, 0.5, 0),
       `Relevant cases` = c("Autosomes in either sex, X chromosomes in males, Z chromosomes in females, chromosomes in haplodiploids of either sex when reproducing sexually, cytoplasmic chromosomes in males",
                            "X chromosomes in females, Z chromosomes in males",
                            "Y chromosomes in males, W chromosomes in females, cytoplasmic chromosomes in females")) %>% 
  pander(split.cell = 20, split.table = Inf)
Prob. that opposite sex gametes carry focal locus Relevant cases
1 Autosomes in either sex, X chromosomes in males, Z chromosomes in females, chromosomes in haplodiploids of either sex when reproducing sexually, cytoplasmic chromosomes in males
0.5 X chromosomes in females, Z chromosomes in males
0 Y chromosomes in males, W chromosomes in females, cytoplasmic chromosomes in females

Once again we can find the condition where breeding with a relative returns greater fitness than an inbreeding avoidance strategy, this time accounting for genetic architecture

\[\delta \lt \frac{ra}{1 + ra}\] Ignoring frequency dependence for now, we can plot the new slopes produced by varying \(r\) and \(a\):

Code
inbreeding_maximum_function_2 <- function(r, a){
  (r*a / (r*a + 1))}

parameters_2 <- expand_grid(r = seq(from = 0, to = 1, by = 0.01),
                          a = c(0, 0.5, 1))

inbreeding_equilibria_2 <- 
  map2_dfr(parameters_2 %>% select(r), 
           parameters_2 %>% select(a), 
           inbreeding_maximum_function_2) %>% 
  rename(depression_threshold = r) %>% 
  bind_cols(parameters_2)

inbreeding_equilibria_2 %>% 
  mutate(a = case_when(a == 0 ~ "a = 0",
                       a == 0.5 ~ "a = 0.5",
                       a == 1 ~ "a = 1")) %>% 
  mutate(a = as.factor(a)) %>% 
  ggplot(aes(x = r, y = depression_threshold, linetype = a, label = a)) +
  geom_textline(linewidth = 0.8, size = 5) + 
 # scale_colour_manual(values = c("0" = met.brewer("Kandinsky", 4)[1], "0.5" =  met.brewer("Kandinsky", 4)[2], "1" = met.brewer("Kandinsky", 4)[3])) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = '_r_, the relatedness coefficient',
       y = ~delta~'(inbreeding depression)',
       linetype = expression(~italic(a)~', the intersex correlation between loci')) +
  scale_x_continuous(expand = c(0.009, 0)) + 
  scale_y_continuous(expand = c(0.025, 0)) +
  theme_bw() +
  theme(text = element_text(size = 18),
        legend.position = "none",
        axis.title.x = element_markdown())

\(~\)

Differences between the sexes beyond genetic architecture

Parker’s seminal equations:

In his 1979 book chapter, Parker considered the inclusive fitness results of breeding with a relative and identified that females and males should have different tolerances for inbreeding depression. The key departure from the unlimited polygyny case presented in the above equations is that a cost to future reproductive success is included for males, i.e. due to finite sperm production, parental care, or harmful mating behaviour such as sexual cannibalism.

For males, Parker found that selection would favour inbreeding with a sister (full-sib) who could otherwise outcross when:

\[n(1 - \delta) + rn(1- \delta) - cn \gt rn\]

the first term is the direct number of alleles propagated, the second term is the indirect number of alleles propagated (note that this is weighted by relatedness), the third term is the direct number of alleles that were not directly propagated by the male through outcrossing, and the final opposing term is the number of alleles that would’ve been transmitted had his sister outcrossed (and he forgone mating).

\(c\) is the cost of the present mating, relative to what is lost for a female. This can be considered a ratio of parental investment. When \(c = 1\) parental investment in the current bout of reproduction is even between the sexes. Alternatively, if males only contribute cheaply produced sperm to an incestuous mating, the cost of mating is likely very small relative to females i.e. \(c ~ 0\).

We add the \(a\) variable to the equation and letting \(n = 1\), simplify to

\[(1-\delta) + ra(1-\delta) - c \gt ra\]

We can again find the condition where breeding with a relative returns greater fitness than an inbreeding avoidance strategy:

\[\delta_{male} = \frac{1 - c}{1 + ra}\]

Parker then modelled the condition for monandrous females to prefer incestuous matings when also presented with an outcrossing opportunity.

\[n(1 - \delta) + rn(1-\delta) - crn \gt n\]

which we can write as

\[(1-\delta) + ra(1-\delta) - rac \gt 1\]

the inbreeding depression threshold is

\[\delta_{female} = \frac{ra - rac}{1 + ra}\]

Note that when \(c = 0\), this is equivalent to the \(\delta\) threshold found in the single mating case.

Plot the relationship between \(r\) and \(\delta\) for several values of \(c\) and \(a\)

Code
Parker_cost_data <- 
  expand_grid(r = seq(from = 0, to = 1, by = 0.01),
              a = c(0, 0.5, 1),
              c = c(0, 0.5, 0.9),
              Sex = c("Female", "Male")) %>% 
  mutate(inbreeding_depression = case_when(Sex == "Female" ~ (r*a - r*a*c) / (r*a + 1),
                                           Sex == "Male" ~ (1 - c) / (r*a + 1)))

  Parker_cost_data %>% 
  mutate(a = as.factor(a)) %>% 
  ggplot(aes(x = r, y = inbreeding_depression, linetype = a, colour = Sex)) +
  geom_line(linewidth = 0.9) + 
  scale_colour_manual(values = c("Female" = met.brewer("Peru1", 6)[2], "Male" =  met.brewer("Peru1", 6)[3])) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(y = expression(delta), 
       x = expression(italic("r")),
       linetype = expression(italic("a"))) +
  scale_x_continuous(expand = c(0.009, 0)) + 
  scale_y_continuous(expand = c(0.01, 0)) +
  facet_wrap(~c, nrow = 3, labeller = label_glue('Male mating investment (c): {c}')) +
  theme_bw() +
  theme(text = element_text(size = 14),
        strip.background = element_rect(fill = "Aliceblue", linewidth = .5))

Plotting sexual and intragenomic conflicts

Code
resolution <- 200

parameters <- 
  expand_grid(
    r = seq(0, 1, length = resolution),
    a = c(0, 0.5, 1),
    c = c(0, 0.25, 0.5, 0.75, 1),
    D = seq(0, 1, length = resolution)) # D represents inbreeding depression)

analytical_results <-
  parameters %>% 
  mutate(female_inbreeding_fitness = (1-D) + (r*a*(1-D)) - (r*a*c),
         female_inbreeding_fitness_cytoplasmic = (1-D) + 0 + 0,
         male_inbreeding_fitness = (1-D) + r*a*(1-D) - c,
         male_inbreeding_fitness_cytoplasmic = 0 + r*1*(1-D) - c,
         female_outbreeding_fitness = 1,
         male_outbreeding_fitness = r*a,
         female_fitness_contrast = female_inbreeding_fitness - female_outbreeding_fitness, # this is close to a selection coefficient
         female_cytoplasmic_fitness_contrast = female_inbreeding_fitness_cytoplasmic - female_outbreeding_fitness,
         male_fitness_contrast = male_inbreeding_fitness - male_outbreeding_fitness,
         male_cytoplasmic_fitness_contrast = male_inbreeding_fitness_cytoplasmic)

When is inbreeding favoured in each sex?

In the figure below, the dashed lines indicate the level of inbreeding depression that can be tolerated for a given value of \(r\). The plot is split into panels by \(c\), the cost of mating for males relative to females and \(a\), the probability that a locus present in one sex is also present in the gametes produced by the other.

Code
pal1 <- met.brewer("OKeeffe1", n=100, direction = -1)
pal2 <- met.brewer("Hiroshige", n=50, direction = -1)
#pal2 <- moma.colors("Avedon", n = 50, direction = 1)

oranges <- c("#ffe6b7", "#ffd06f", "#f7aa58", "#ef8a47", "#e76254") 

blues <- c("#aadce0", "#72bcd5", "#528fad", "#376795", "#1e466e")

# Lowest colour for all gradients
gradient_base <- oranges[1]

gradient_max <- oranges[5]

# Create a list of gradients for each colour 2 to 10 over five steps from 
# gradient_base grey (low) to colour (high)
my_gradient <- colorRampPalette(c(gradient_base, gradient_max))(50)
  
Female_plot <-
  analytical_results %>%
  filter(c == 0 | c == 0.5) %>% 
  ggplot(aes(x = r, y = D)) +
  geom_blank() +
  geom_raster(aes(fill = female_fitness_contrast)) + 
  stat_contour(aes(z = female_fitness_contrast), colour = "black", binwidth = 0.25,
               breaks = c(-1.25, -1, -0.75, -0.5, -0.25, 0.25, 0.5, 0.75, 1, 1.25)) +
  stat_contour(aes(z = female_fitness_contrast), colour = "black", breaks = 0,
               linetype = 2) +
  scale_fill_gradientn(colours = pal1, breaks = c(-2, -1, 0, 1, 2), limits = c(-2, 2)) +
  facet_wrap(c ~ a, 
             scales = "free", nrow = 2, strip.position = c("top"),
             labeller = label_glue('Prob. related gamete carries inbreeding locus = {`a`}\nMale mating cost = {`c`}')) +
  labs(x = '_r_, the relatedness coefficient',
       y = ~delta~'(inbreeding depression)',
       fill = "Inbreeding fitness",
       title = "A. Alleles present in females") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 7),
        axis.title.x = element_markdown(),
        plot.title = element_text(hjust = 0.5))


Male_plot <-
  analytical_results %>%
  filter(c == 0 | c == 0.5) %>% 
  ggplot(aes(x = r, y = D)) +
  geom_blank() +
  geom_raster(aes(fill = male_fitness_contrast)) + 
  stat_contour(aes(z = male_fitness_contrast), colour = "black", binwidth = 0.25,
              breaks = c(-1.25, -1, -0.75, -0.5, -0.25, 0.25, 0.5, 0.75, 1, 1.25)) +
  stat_contour(aes(z = male_fitness_contrast), colour = "black", breaks = 0,
               linetype = 2) +
   scale_fill_gradientn(colours = pal1, breaks = c(-2, -1, 0, 1, 2), limits = c(-2, 2)) +
   facet_wrap(c ~ a, 
             scales = "free", nrow = 2, strip.position = c("top"),
             labeller = label_glue('Prob. related gamete carries inbreeding locus = {`a`}\nMale mating cost = {`c`}')) +
  labs(x = '_r_, the relatedness coefficient',
       y = ~delta~'(inbreeding depression)',
       fill = "Inbreeding fitness",
       title = "B. Alleles present in males") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", linewidth = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 7),
        axis.title.x = element_markdown(),
        plot.title = element_text(hjust = 0.5))

Female_plot / Male_plot + plot_layout(guides = "collect")

When is there evolutionary conflict over inbreeding?

To quantify conflict in its various forms we use a version of the I index presented in Innocenti and Morrow (2010).

Code
Intragenomic_conflict_females <-
  analytical_results %>%
  select(1:4, female_fitness_contrast) %>%
  pivot_wider(names_from = a, values_from = female_fitness_contrast) %>% 
  #mutate(`Autosome - X` = `1` - `0.5`,
   #      `Autosome - Z` = 0,
    #     `Autosome - W & Z - W` = `1` - `0`) %>% 
  mutate(`Autosome vs X` = `1` * `0.5` / sqrt(((`1`)^2 + (`0.5`)^2)/2),
         `Autosome vs Z` = `1` * `1` / sqrt(((`1`)^2 + (`1`)^2)/2),
         `Autosome & Z vs W` = `1` * `0` / sqrt(((`1`)^2 + (`0`)^2)/2)) %>% 
  pivot_longer(cols = 7:9, names_to = "contrast", values_to = "Evolutionary_conflict") %>% 
  mutate(relationship = case_when(
    contrast == "Autosome vs X" & `1` > 0 & `0.5` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "Autosome vs X" & `1` < 0 & `0.5` < 0 ~
       "Inbreeding deleterious in both contexts",
    
    contrast == "Autosome vs Z" & `1` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "Autosome vs Z" & `1` < 0 ~
       "Inbreeding deleterious in both contexts",
    
    contrast == "Autosome & Z vs W" & `1` > 0 & `0` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "Autosome & Z vs W" & `1` < 0 & `0` < 0 ~
       "Inbreeding deleterious in both contexts",
    `0` == 0 & `0.5` == 0 & `1` == 0 ~ "Inbreeding deleterious in both contexts",  
    .default = "Intragenomic conflict")) %>% 
  mutate(contrast = fct_relevel(contrast, "Autosome vs X", "Autosome vs Z",
                                "Autosome & Z vs W")) %>% 
  filter(c == 0)  # remove if we want more c values

Intragenomic_conflict_males <-
  analytical_results %>%
  select(1:4, male_fitness_contrast) %>%
  pivot_wider(names_from = a, values_from = male_fitness_contrast) %>% 
  #mutate(`Autosome - X` = 0,
   #      `Autosome - Z` = `1` - `0.5`,
    #     `Autosome - Y & X - Y` = `1` - `0`) %>%
  mutate(`X vs Autosome` = `1` * `1` / sqrt(((`1`)^2 + (`1`)^2)/2),
         `Z vs Autosome` = `1` * `0.5` / sqrt(((`1`)^2 + (`0.5`)^2)/2),
         `Y vs Autosome & X` = `1` * `0` / sqrt(((`1`)^2 + (`0`)^2)/2)) %>% 
  pivot_longer(cols = 7:9, names_to = "contrast", values_to = "Evolutionary_conflict") %>% 
  mutate(relationship = case_when(
    contrast == "X vs Autosome" & `1` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "X vs Autosome" & `1` < 0 ~
      "Inbreeding deleterious in both contexts",
    
    contrast == "Z vs Autosome" & `1` > 0 & `0.5` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "Z vs Autosome" & `1` < 0 & `0.5` < 0 ~
      "Inbreeding deleterious in both contexts",
    
    contrast == "Y vs Autosome & X" & `1` > 0 & `0` > 0 ~
      "Inbreeding favoured in both contexts",
    contrast == "Y vs Autosome & X" & `1` < 0 & `0` < 0 ~
      "Inbreeding deleterious in both contexts",
    `0` == 0 & `0.5` == 0 & `1` == 0 ~ "Inbreeding deleterious in both contexts",  
    .default = "Intragenomic conflict")) %>% 
  mutate(contrast = fct_relevel(contrast, "X vs Autosome", "Z vs Autosome",
                                "Y vs Autosome & X")) %>% 
  filter(c == 0) # remove if we want more c values
    

make_genomic_conflict_plot <- 
  function(data, enter_title, colour_pal){
    data %>% 
      ggplot(aes(x = r, y = D)) +
      geom_blank() +
      geom_tile(data = data %>% filter(relationship == "Intragenomic conflict"),
                aes(fill = Evolutionary_conflict*-1)) + 
      #geom_tile(data = data,
      #         aes(fill = Evolutionary_conflict)) + 
      scale_fill_gradientn(colours = colour_pal, limits = c(0, 0.6), #na.value = "white",
                           breaks = c(0, 0.6),
                           labels = c("Weaker conflict", "Stronger conflict")) +
      labs(fill = "Evolutionary conflict") +
      new_scale_fill() +
      geom_tile(data = data %>% filter(relationship != "Intragenomic conflict"),
                aes(fill = relationship), alpha = 0.75) +
      scale_fill_manual(values = c("#fbe6c5", "#d2fbd4"), 
                        labels = c("Inbreeding deleterious\nin both contexts", 
                               "Inbreeding favoured\nin both contexts")) +
      stat_contour(aes(z = Evolutionary_conflict*-1), colour = "black", binwidth = 25,
                   breaks = c(0, 0.2, 0.4, 0.6)) +
      facet_wrap(~contrast, nrow = 3,
                 scales = "free", strip.position = c("top"),
                 labeller = label_glue('{`contrast`}')) +
      labs(x = '_r_ (relatedness coefficient)',
           y = ~delta~'(inbreeding depression)',
           fill = "Evolutionary concordance",
           title = enter_title) +
      scale_x_continuous(expand = c(0, 0)) + 
      scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
      theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
            panel.grid.minor = element_blank(),
            strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
            axis.title.x = element_markdown(),
            plot.title = element_text(hjust = 0.5, size = 7))
  }

icf <- make_genomic_conflict_plot(Intragenomic_conflict_females, "A. Female inter-chromosomal conflict", oranges)
icm <- make_genomic_conflict_plot(Intragenomic_conflict_males, "B. Male inter-chromosomal conflict", oranges)

Make column 3

Code
autosomal_data <-
  analytical_results %>%
  filter(a == 1, c != 1, c != 0.25) %>% 
  #mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>%
  mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast / 
           sqrt(((female_fitness_contrast)^2 + (male_fitness_contrast)^2)/2)) %>% 
  mutate(relationship = 
           case_when(female_fitness_contrast < 0 & male_fitness_contrast > 0 
                     ~ "Sexual conflict",
                     female_fitness_contrast < 0 & male_fitness_contrast < 0 
                     ~ "Inbreeding deleterious in both contexts",
                     female_fitness_contrast > 0 & male_fitness_contrast > 0 
                     ~ "Inbreeding favoured in both contexts"),
         Location = "Autosome") 

X_data <- 
  analytical_results %>% 
  filter(a == 0.5, c != 1, c != 0.25) %>%  
  select(1:4, contains("female")) %>% 
  rename(a_female = a) %>% # this step makes the join work as intended
  left_join(
    analytical_results %>% 
      filter(a == 1, c != 1, c != 0.25) %>% 
      select(1:4, starts_with("male")) %>% 
      rename(a_male = a) # this step makes the join work as intended
  ) %>%
  #mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>% 
  mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast / 
           sqrt(((female_fitness_contrast)^2 + (male_fitness_contrast)^2)/2)) %>% 
  mutate(relationship = 
           case_when(female_fitness_contrast < 0 & male_fitness_contrast > 0 
                     ~ "Sexual conflict",
                     female_fitness_contrast < 0 & male_fitness_contrast < 0 
                     ~ "Inbreeding deleterious in both contexts",
                     female_fitness_contrast > 0 & male_fitness_contrast > 0 
                     ~ "Inbreeding favoured in both contexts"),
         Location = "X")

Z_data <-
  analytical_results %>% 
  filter(a == 1, c != 1, c != 0.25) %>%  
  select(1:4, contains("female")) %>% 
  rename(a_female = a) %>% 
  left_join(
    analytical_results %>% 
      filter(a == 0.5, c != 1, c != 0.25) %>% 
      select(1:4, starts_with("male")) %>% 
      rename(a_male = a)
  ) %>% 
  #mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>%
  mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast / 
           sqrt(((female_fitness_contrast)^2 + (male_fitness_contrast)^2)/2)) %>% 
  mutate(relationship = 
           case_when(female_fitness_contrast < 0 & male_fitness_contrast > 0 
                     ~ "Sexual conflict",
                     female_fitness_contrast < 0 & male_fitness_contrast < 0 
                     ~ "Inbreeding deleterious in both contexts",
                     female_fitness_contrast > 0 & male_fitness_contrast > 0 
                     ~ "Inbreeding favoured in both contexts"),
         Location = "Z")

plotting_data <- bind_rows(autosomal_data, X_data, Z_data) %>% filter(c == 0)

Sexual_conflict_plot <-
  plotting_data %>%
  ggplot(aes(x = r, y = D)) +
  geom_blank() +
  geom_tile(data = plotting_data %>% filter(relationship == "Sexual conflict"),
            aes(fill = Evolutionary_conflict*-1)) + 
  #scale_fill_gradientn(colours = pal2, limits = c(-1.2, 1.1), #na.value = "white",
  #                    labels = c("Strong conflict, female (+)", -0.5, 
  #                              "No conflict", 0.5, "Strong conflict, male (+)")) +
  scale_fill_gradientn(colours = oranges, limits = c(0, 0.6), #na.value = "white",
                       breaks = c(0, 0.6),
                       labels = c("Weaker conflict", "Stronger conflict")) +
  labs(fill = "Evolutionary conflict") +
  new_scale_fill() +
  geom_tile(data = plotting_data %>% filter(relationship != "Sexual conflict"),
            aes(fill = relationship), alpha = 0.75) +
  scale_fill_manual(values = c("#fbe6c5", "#d2fbd4"), 
                    labels = c("Inbreeding deleterious\nin both contexts", 
                               "Inbreeding favoured\nin both contexts")) +
  stat_contour(aes(z = Evolutionary_conflict*-1), colour = "black", binwidth = 25,
               breaks = c(0, 0.2, 0.4, 0.6)) +
  facet_wrap(~Location, nrow = 3,
             scales = "free", strip.position = c("top"),
             labeller = label_glue('{`Location`}')) +
  labs(x = '_r_ (relatedness coefficient)',
       y = ~delta~'(inbreeding depression)',
       fill = "Evolutionary concordance",
       title = "C. Intra-chromosomal sexual conflict") +
  scale_x_continuous(expand = c(0, 0)) + 
  scale_y_continuous(expand = c(0, 0)) + # labels = c(0, 25, 50, 75, 90)) +
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.title.x = element_markdown(),
        plot.title = element_text(hjust = 0.5, size = 7))

Make Figure 1

Code
icf + icm + Sexual_conflict_plot + plot_layout(guides = "collect")

Panels in columns A and B show the parameter space where inbreeding increases the propagation of one type of chromosome, but has an opposite, deleterious affect on the propagation of a second chromosome class - so called inter-chromosomal conflict.

Column A shows regions of inter-chromosomal conflict within females, where chromosome classes with a high probability of encountering gametes carrying the inbreeding locus have high fitness in the conflict zone. In females, autosomal and Z chromosomes are aligned in the fitness outcomes of inbreeding; these cases reflect the implicit assumptions of previous models exploring the inclusive fitness effects of inbreeding. In some cases, the fitness interests of chromosome classes align e.g. autosomes and Z chromosomes (this is only true for females).

The second column shows inter-chromosomal in males. Importantly, the direction of selection for inbreeding is reversed in comparison to the female case in the zones of conflict. In males, chromosome classes with high probabilities of inbreeding loci being found in female gametes now have low fitness in the conflict space.

Panels in column C show where loci are expected to be under intra-chromosomal sexual conflict over inbreeding preference. In regions of sexual conflict, inbreeding preference is always favoured in males, but has negative fitness consequences if expressed by females (assuming that males invest less into mating than females). Note that intra-chromosomal conflict encompasses both the intra- and inter-locus forms of sexual conflict.

Make Figure 2

We can also plot the cost to female direct reproductive output,

Code
fitness_loss_data <-
  tibble(chromosome = c("Autosomal", "X", "Z", "Y", "W", "Cytoplasmic"),
       a_female = c(1, 0.5, 1, NA, 0, 0),
       a_male = c(1, 1, 0.5, 0, NA, 1)) %>% 
  mutate("delta[female]~when~r==0.5" = 0.5*a_female / (1 + 0.5*a_female),
         "delta[male]~when~r==0.5" = 1 / (1 + 0.5*a_male),
         "delta[male]~when~r==0.5" = case_when(chromosome == "Cytoplasmic" ~ 0,
                                               .default = `delta[male]~when~r==0.5`)) %>%
  pivot_longer(cols = 4:5, names_to = "Sex", values_to = "Inbreeding tolerance") %>% 
  filter(!is.na(`Inbreeding tolerance`)) 

fitness_loss_f_delta <-
  fitness_loss_data %>% 
  filter(Sex == "delta[female]~when~r==0.5") %>% 
  mutate(chromosome = fct_relevel(chromosome, c("Autosomal", "X", "Z", "W", "Cytoplasmic"))) %>% 
  ggplot(aes(x = chromosome, y = 1 - `Inbreeding tolerance`)) +
  geom_col(fill = "#d0e2af", alpha = 0.75) +
  facet_wrap(~Sex, labeller = label_parsed) +
  scale_y_continuous(expand = c(0.0, 0.0), limits = c(0, 1.01), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  labs(x = "Genomic location",
       y = "Female organismal-level direct fitness") +
  theme_bw() +
  theme(strip.background = element_rect(colour = "black", fill = "#d0e2af"),
        strip.text = element_text(size = 12),
        axis.title = element_text(size = 12))

fitness_loss_m_delta <-
  fitness_loss_data %>% 
  filter(Sex == "delta[male]~when~r==0.5") %>% 
  mutate(chromosome = fct_relevel(chromosome, c("Autosomal", "X", "Z", "Y", "Cytoplasmic"))) %>% 
  ggplot(aes(x = chromosome, y = 1 - `Inbreeding tolerance`)) +
  geom_col(fill = "#d0e2af", alpha = 0.75) +
    facet_wrap(~Sex, labeller = label_parsed) +
  scale_y_continuous(expand = c(0.0, 0.0), limits = c(0, 1.01), breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) + 
  labs(x = "Genomic location",
       y = "Female organismal-level direct fitness") +
  theme_bw() +
  theme(strip.background = element_rect(colour = "black", fill = "#d0e2af"),
        strip.text = element_text(size = 12),
        axis.title = element_text(size = 12))

fitness_loss_f_delta + fitness_loss_m_delta +
  plot_layout(axis_titles = "collect")

Figure 2. direct female reproductive output at the threshold where chromosome-specific selection should no longer favour inbreeding. Reproductive output is expressed relative to females that exclusively outcross.

The simulation

Parker’s \(c\) value offers a simple, intuitive way to model the cost of mating for males relative to females. However, other methods better capture the dynamics of real populations, where both sexes also run the risk of going unmated. To incorporate both costs of mating and matelessness, we simulate the invasion of an allele that encodes a preference for inbreeding, for loci on various chromosomes that are expressed in either females or males.

The simulation progresses through time by jumping from event to event. Events are calculated via a gillespie-like algorithm, which involves both memory-less processes and others that are maintained through time. Events can trigger state changes for the individuals in the population, leading to death, mating, inheritance of a breeding site or offspring production. The simulation ends when the inbreeding allele is purged from the population, it reaches a high frequency in the population, the population goes extinct, or the time limit elapses.

Following Ekrem and Kokko (2023), we find the lifespans of N diploid individuals in a sexually reproducing population to initialise the simulation. Mortality events are drawn from an exponential distribution with rate \(\lambda = 1/\mu\), where \(\mu\) is the mean time till death. The probability of mortality is therefore constant across life for each sex.

Lets have a look at the mean lifespan drawn from the exponential distribution parameterised with different mortality rates:

Code
#exponential_draws <-
  tibble(`Mean lifespan` = seq(from = 0.001, to = 1.001, by = 0.1)) %>% 
  mutate(`Corresponding mortality rate` = 1/`Mean lifespan`) %>% 
  arrange(-`Mean lifespan`) %>% 
  pander()
Mean lifespan Corresponding mortality rate
1.001 0.999
0.901 1.11
0.801 1.248
0.701 1.427
0.601 1.664
0.501 1.996
0.401 2.494
0.301 3.322
0.201 4.975
0.101 9.901
0.001 1000

We use these values to inform our parameterisation of inbreeding depression (\(\delta\)). Individuals that don’t inbreed produce offspring with a mean lifespan of 1 time unit (\(\lambda = 1\)). For inbred offspring, we model inbreeding depression by increasing the risk of mortality (\(\lambda > 1\)), lowering the mean lifespan of inbred individuals relative to their outbred counterparts.

Each individual possesses an inbreeding locus with two possible alleles. The A_O allele encodes inbreeding avoidance, whereas the A_I allele encodes inbreeding preference. The number_mutants argument allows us to set the number of mutations that the population is originally seeded with.

Female-male encounters are found using the mass action law (also called fertilisation kinetics in this biological context), with male search efficiency \(v\). To get an idea of how \(v\) influences mating encounters, we plot the population wide male-female encounter timestamps. The dashed line shows the mean lifespan of outbred individuals in our simulations. Note that these are just encounters; whether they lead to mating depends upon past mating experiences of both females and males.

Code
data <-
  expand_grid(v = c(0.05, 0.5, 2, 5),
              lifespan = c(0.01, 0.1, 0.5, 1),
         pop_size = rep(100, 10000)) %>% 
  mutate(encounter = rexp(160000, rate = v)) %>% 
  group_by(v, lifespan) %>% 
  summarise(n_females = n(),
              lifetime_encounters = sum(encounter < lifespan),
              prop_females_encountered = lifetime_encounters / n_females) %>% 
  ungroup() %>% 
   mutate(v = as.factor(v),
         lifespan = as.factor(lifespan))

data  %>% 
  ggplot(aes(x = lifespan, y = prop_females_encountered)) +
   geom_col(fill = "#d0e2af", alpha = 1) +
  scale_y_continuous(expand = c(0.0, 0.01), limits = c(0, 1)) + 
  labs(x = "Mean lifespan",
       y = "Proportion of female population encountered throughout lifetime") +
  #facet_wrap(~ lifespan) +
  facet_wrap(~v, nrow = 2,
             scales = "free", strip.position = c("top"),
             labeller = label_glue('Male search efficiency = {`v`}')) +
  theme_bw() +
  theme(strip.background = element_rect(colour = "black", fill = "aliceblue"),
        strip.text = element_text(size = 12),
        axis.title = element_text(size = 12))

Figure. SX: the proportion of the female population encountered by a single male across his lifetime, split by how efficiently males search for females. If there are 100 females in the population, a value of 0.1 indicates the male encounters 10 females on average. The decline in potential mating encounters as lifespan decreases indicates the severity of the inbreeding depression cost for males.

The cost of inbreeding depression manifests differently for females, where the ability to reproduce depends on mating with a male and securing a breeding site. If the population is at carrying capacity, the latter criterion is a function of age, as breeding sites are allocated randomly to non-breeding females once a breeder dies. Therefore, the longer a female lives, the more chances she has of becoming a breeder. If the population is below carrying capacity, the cost of inbreeding is less severe, as a new female can immediately secure a breeding site, leaving mating as the only restriction to offspring production.

Build a mating table

To simulate the invasion of alleles found on different chromosomes, we consider simple single-locus, two-allele dynamics. For each chromosomal inheritance system, we build a mating table, which tracks the production of new genotypes from each possible mother-father pair of genotypes in the population. Inheritance is mendelian, except for cytoplamsic chromosomes, which are exclusively maternally inherited.

The simulation tracks the invasion of an allele encoding a preference for inbreeding (denoted I), into a population initially dominated by an allele that encodes inbreeding avoidance (denoted O for outbreeding).

Code
make_mating_table <- function(gene_location){
  
  make_offspring <- function(X, Y, type, zygote_freq, gene_location){
    tibble(Female_genotype = X,
           Male_genotype = Y,
           type,
           zygote_freq,
           locus_type = gene_location)
  }
  
  # Specify the possible offspring genotypes for all the potential crosses; we use these for the type argument in the make_offspring function
  
  # autosomal
  
  # II x II
  
  a_genotype_1 <- c("A_IA_I.Female", "A_IA_I.Male")
  
  # II x IO
  
  a_genotype_2 <- c("A_IA_I.Female", "A_IA_I.Male", 
                    "A_IA_O.Female", "A_IA_O.Male")
  
  # II x OO
  
  a_genotype_3 <- c("A_IA_O.Female", "A_IA_O.Male")
  
  # IO x IO
  
  a_genotype_4 <- c("A_IA_I.Female", "A_IA_I.Male", 
                    "A_IA_O.Female", "A_IA_O.Male", 
                    "A_OA_O.Female", "A_OA_O.Male")
  
  # IO x OO
  
  a_genotype_5 <- c("A_IA_O.Female", "A_IA_O.Male",
                    "A_OA_O.Female", "A_OA_O.Male")
  
  # OO x OO
  
  a_genotype_6 <- c("A_OA_O.Female", "A_OA_O.Male")
  
  # XY
  
  
  # II x IY_I
  
  xy_genotype_1 <- c("X_IX_I.Female", "X_IY_I.Male")
  
  # II x IY_O
  
  xy_genotype_2 <- c("X_IX_I.Female", "X_IY_O.Male")
  
  # II x OY_I
  
  xy_genotype_3 <- c("X_IX_O.Female", "X_IY_I.Male")
  
  # II x OY_O
  
  xy_genotype_4 <- c("X_IX_O.Female", "X_IY_O.Male")
  
  # IO x IY_I
  
  xy_genotype_5 <- c("X_IX_I.Female", "X_IY_I.Male",
                     "X_IX_O.Female", "X_OY_I.Male")
  
  # IO x IY_O
  
  xy_genotype_6 <- c("X_IX_I.Female", "X_IY_O.Male", 
                     "X_IX_O.Female", "X_OY_O.Male")
  
  # IO x OY_I
  
  xy_genotype_7 <- c("X_IX_O.Female", "X_IY_I.Male",
                     "X_OX_O.Female", "X_OY_I.Male")
  
  # IO x OY_O
  
  xy_genotype_8 <- c("X_IX_O.Female", "X_IY_O.Male",
                     "X_OX_O.Female", "X_OY_O.Male")
  
  # OO x IY_I
  
  xy_genotype_9 <- c("X_IX_O.Female", "X_OY_I.Male")
  
  # OO x IY_O
  
  xy_genotype_10 <- c("X_IX_O.Female", "X_OY_O.Male")
  
  # OO x OY_I
  
  xy_genotype_11 <- c("X_OX_O.Female", "X_OY_I.Male")
  
  # OO x OY_O
  
  xy_genotype_12 <- c("X_OX_O.Female", "X_OY_O.Male")
  
  # ZW
  
  # IW_I x II
  
  zw_genotype_1 <- c("Z_IZ_I.Male", "Z_IW_I.Female")
  
  # IW_I x IO
  
  zw_genotype_2 <- c("Z_IZ_I.Male", "Z_IZ_O.Male", 
                     "Z_IW_I.Female", "Z_OW_I.Male")
  
  # IW_I x OO
  
  zw_genotype_3 <- c("Z_IZ_O.Male", "Z_OW_I.Female")
  
  # IW_O x II
  
  zw_genotype_4 <- c("Z_IZ_I.Male", "Z_IW_O.Female")
  
  # IW_O x IO
  
  zw_genotype_5 <- c("Z_IZ_I.Male", "Z_IZ_O.Male",
                     "Z_IW_O.Female", "Z_OW_O.Female")
  
  # IW_O x OO
  
  zw_genotype_6 <- c("Z_IZ_O.Male", "Z_OW_O.Female")
  
  # OW_I X II
  
  zw_genotype_7 <- c("Z_IZ_O.Male", "Z_IW_I.Female")
  
  # OW_I x IO
  
  zw_genotype_8 <- c("Z_IZ_O.Male", "Z_OZ_O.Male",
                     "Z_IW_I.Female", "Z_OW_I.Female")
  
  # OW_I x OO
  
  zw_genotype_9 <- c("Z_OZ_O.Male", "Z_OW_I.Female")
  
  # OW_O X II
  
  zw_genotype_10 <- c("Z_IZ_O.Male", "Z_IW_O.Female")
  
  # OW_O x IO
  
  zw_genotype_11 <- c("Z_IZ_O.Male", "Z_OZ_O.Male",
                      "Z_IW_O.Female", "Z_OW_O.Female")
  
  # OW_O x OO
  
  zw_genotype_12 <- c("Z_OZ_O.Male", "Z_OW_O.Female")
  
  # cytoplasmic
  
  # I x I
  # I x O
  
  c_genotype_1 <- c("C_I.Female", "C_I.Male")
  
  # O x O
  # O x I
  
  c_genotype_2 <- c("C_O.Female", "C_O.Male")
  
  
  
  # Now calculate the zygote frequencies for each cross
  
  # autosomal
  
  # even frequency of two offspring genotypes
  
  freq_2 <- rep(0.5, 2)
  
  # even frequency between four offspring types
  
  freq_4 <- rep(0.25, 4)
  
  # when there are 6 offspring genotypes
  
  freq_6 <- c(0.125, 0.125,
              0.25, 0.25,
              0.125, 0.125)
  
  bind_rows(
    list(
      make_offspring("A_IA_I", "A_IA_I", a_genotype_1, freq_2, "autosomal"),
      make_offspring("A_IA_I", "A_IA_O", a_genotype_2, freq_4, "autosomal"),
      make_offspring("A_IA_I", "A_OA_O", a_genotype_3, freq_2, "autosomal"),
      make_offspring("A_IA_O", "A_IA_I", a_genotype_2, freq_4, "autosomal"),
      make_offspring("A_IA_O", "A_IA_O", a_genotype_4, freq_6, "autosomal"),
      make_offspring("A_IA_O", "A_OA_O", a_genotype_5, freq_4, "autosomal"),
      make_offspring("A_OA_O", "A_IA_I", a_genotype_3, freq_2, "autosomal"),
      make_offspring("A_OA_O", "A_IA_O", a_genotype_5, freq_4, "autosomal"),
      make_offspring("A_OA_O", "A_OA_O", a_genotype_6, freq_2, "autosomal"),
      
      make_offspring("X_IX_I", "X_IY_I", xy_genotype_1, freq_2, "XY"),
      make_offspring("X_IX_I", "X_IY_O", xy_genotype_2, freq_2, "XY"),
      make_offspring("X_IX_I", "X_OY_I", xy_genotype_3, freq_2, "XY"),
      make_offspring("X_IX_I", "X_OY_O", xy_genotype_4, freq_2, "XY"),
      make_offspring("X_IX_O", "X_IY_I", xy_genotype_5, freq_4, "XY"),
      make_offspring("X_IX_O", "X_IY_O", xy_genotype_6, freq_4, "XY"),
      make_offspring("X_IX_O", "X_OY_I", xy_genotype_7, freq_4, "XY"),
      make_offspring("X_IX_O", "X_OY_O", xy_genotype_8, freq_4, "XY"),
      make_offspring("X_OX_O", "X_IY_I", xy_genotype_9, freq_2, "XY"),
      make_offspring("X_OX_O", "X_IY_O", xy_genotype_10, freq_2, "XY"),
      make_offspring("X_OX_O", "X_OY_I", xy_genotype_11, freq_2, "XY"),
      make_offspring("X_OX_O", "X_OY_O", xy_genotype_12, freq_2, "XY"),
      
      make_offspring("Z_IW_I", "Z_IZ_I", zw_genotype_1, freq_2, "ZW"),
      make_offspring("Z_IW_I", "Z_IZ_O", zw_genotype_2, freq_4, "ZW"),
      make_offspring("Z_IW_I", "Z_OZ_O", zw_genotype_3, freq_2, "ZW"),
      make_offspring("Z_IW_O", "Z_IZ_I", zw_genotype_4, freq_2, "ZW"),
      make_offspring("Z_IW_O", "Z_IZ_O", zw_genotype_5, freq_4, "ZW"),
      make_offspring("Z_IW_O", "Z_OZ_O", zw_genotype_6, freq_2, "ZW"),
      make_offspring("Z_OW_I", "Z_IZ_I", zw_genotype_7, freq_2, "ZW"),
      make_offspring("Z_OW_I", "Z_IZ_O", zw_genotype_8, freq_4, "ZW"),
      make_offspring("Z_OW_I", "Z_OZ_O", zw_genotype_9, freq_2, "ZW"),
      make_offspring("Z_OW_O", "Z_IZ_I", zw_genotype_10, freq_2, "ZW"),
      make_offspring("Z_OW_O", "Z_IZ_O", zw_genotype_11, freq_4, "ZW"),
      make_offspring("Z_OW_O", "Z_OZ_O", zw_genotype_12, freq_2, "ZW"),
      
      make_offspring("C_I", "C_I", c_genotype_1, freq_2, "cytoplasmic"),
      make_offspring("C_I", "C_O", c_genotype_1, freq_2, "cytoplasmic"),
      make_offspring("C_O", "C_I", c_genotype_2, freq_2, "cytoplasmic"),
      make_offspring("C_O", "C_O", c_genotype_2, freq_2, "cytoplasmic")
    )) %>% 
    filter(locus_type == gene_location)
}
Code
make_mating_table(gene_location = "autosomal") %>% 
  select(locus_type, everything()) %>% 
  rename(zygote_type = type) %>% 
  kable() %>% 
  kable_styling() %>% 
  scroll_box(height = "500px")
locus_type Female_genotype Male_genotype zygote_type zygote_freq
autosomal A_IA_I A_IA_I A_IA_I.Female 0.500
autosomal A_IA_I A_IA_I A_IA_I.Male 0.500
autosomal A_IA_I A_IA_O A_IA_I.Female 0.250
autosomal A_IA_I A_IA_O A_IA_I.Male 0.250
autosomal A_IA_I A_IA_O A_IA_O.Female 0.250
autosomal A_IA_I A_IA_O A_IA_O.Male 0.250
autosomal A_IA_I A_OA_O A_IA_O.Female 0.500
autosomal A_IA_I A_OA_O A_IA_O.Male 0.500
autosomal A_IA_O A_IA_I A_IA_I.Female 0.250
autosomal A_IA_O A_IA_I A_IA_I.Male 0.250
autosomal A_IA_O A_IA_I A_IA_O.Female 0.250
autosomal A_IA_O A_IA_I A_IA_O.Male 0.250
autosomal A_IA_O A_IA_O A_IA_I.Female 0.125
autosomal A_IA_O A_IA_O A_IA_I.Male 0.125
autosomal A_IA_O A_IA_O A_IA_O.Female 0.250
autosomal A_IA_O A_IA_O A_IA_O.Male 0.250
autosomal A_IA_O A_IA_O A_OA_O.Female 0.125
autosomal A_IA_O A_IA_O A_OA_O.Male 0.125
autosomal A_IA_O A_OA_O A_IA_O.Female 0.250
autosomal A_IA_O A_OA_O A_IA_O.Male 0.250
autosomal A_IA_O A_OA_O A_OA_O.Female 0.250
autosomal A_IA_O A_OA_O A_OA_O.Male 0.250
autosomal A_OA_O A_IA_I A_IA_O.Female 0.500
autosomal A_OA_O A_IA_I A_IA_O.Male 0.500
autosomal A_OA_O A_IA_O A_IA_O.Female 0.250
autosomal A_OA_O A_IA_O A_IA_O.Male 0.250
autosomal A_OA_O A_IA_O A_OA_O.Female 0.250
autosomal A_OA_O A_IA_O A_OA_O.Male 0.250
autosomal A_OA_O A_OA_O A_OA_O.Female 0.500
autosomal A_OA_O A_OA_O A_OA_O.Male 0.500
Code
make_mating_table(gene_location = "XY") %>% 
  select(locus_type, everything()) %>% 
  rename(zygote_type = type) %>% 
  kable() %>% 
  kable_styling() %>% 
  scroll_box(height = "500px")
locus_type Female_genotype Male_genotype zygote_type zygote_freq
XY X_IX_I X_IY_I X_IX_I.Female 0.50
XY X_IX_I X_IY_I X_IY_I.Male 0.50
XY X_IX_I X_IY_O X_IX_I.Female 0.50
XY X_IX_I X_IY_O X_IY_O.Male 0.50
XY X_IX_I X_OY_I X_IX_O.Female 0.50
XY X_IX_I X_OY_I X_IY_I.Male 0.50
XY X_IX_I X_OY_O X_IX_O.Female 0.50
XY X_IX_I X_OY_O X_IY_O.Male 0.50
XY X_IX_O X_IY_I X_IX_I.Female 0.25
XY X_IX_O X_IY_I X_IY_I.Male 0.25
XY X_IX_O X_IY_I X_IX_O.Female 0.25
XY X_IX_O X_IY_I X_OY_I.Male 0.25
XY X_IX_O X_IY_O X_IX_I.Female 0.25
XY X_IX_O X_IY_O X_IY_O.Male 0.25
XY X_IX_O X_IY_O X_IX_O.Female 0.25
XY X_IX_O X_IY_O X_OY_O.Male 0.25
XY X_IX_O X_OY_I X_IX_O.Female 0.25
XY X_IX_O X_OY_I X_IY_I.Male 0.25
XY X_IX_O X_OY_I X_OX_O.Female 0.25
XY X_IX_O X_OY_I X_OY_I.Male 0.25
XY X_IX_O X_OY_O X_IX_O.Female 0.25
XY X_IX_O X_OY_O X_IY_O.Male 0.25
XY X_IX_O X_OY_O X_OX_O.Female 0.25
XY X_IX_O X_OY_O X_OY_O.Male 0.25
XY X_OX_O X_IY_I X_IX_O.Female 0.50
XY X_OX_O X_IY_I X_OY_I.Male 0.50
XY X_OX_O X_IY_O X_IX_O.Female 0.50
XY X_OX_O X_IY_O X_OY_O.Male 0.50
XY X_OX_O X_OY_I X_OX_O.Female 0.50
XY X_OX_O X_OY_I X_OY_I.Male 0.50
XY X_OX_O X_OY_O X_OX_O.Female 0.50
XY X_OX_O X_OY_O X_OY_O.Male 0.50
Code
make_mating_table(gene_location = "ZW") %>% 
  select(locus_type, everything()) %>% 
  rename(zygote_type = type) %>% 
  kable() %>% 
  kable_styling() %>% 
  scroll_box(height = "500px")
locus_type Female_genotype Male_genotype zygote_type zygote_freq
ZW Z_IW_I Z_IZ_I Z_IZ_I.Male 0.50
ZW Z_IW_I Z_IZ_I Z_IW_I.Female 0.50
ZW Z_IW_I Z_IZ_O Z_IZ_I.Male 0.25
ZW Z_IW_I Z_IZ_O Z_IZ_O.Male 0.25
ZW Z_IW_I Z_IZ_O Z_IW_I.Female 0.25
ZW Z_IW_I Z_IZ_O Z_OW_I.Male 0.25
ZW Z_IW_I Z_OZ_O Z_IZ_O.Male 0.50
ZW Z_IW_I Z_OZ_O Z_OW_I.Female 0.50
ZW Z_IW_O Z_IZ_I Z_IZ_I.Male 0.50
ZW Z_IW_O Z_IZ_I Z_IW_O.Female 0.50
ZW Z_IW_O Z_IZ_O Z_IZ_I.Male 0.25
ZW Z_IW_O Z_IZ_O Z_IZ_O.Male 0.25
ZW Z_IW_O Z_IZ_O Z_IW_O.Female 0.25
ZW Z_IW_O Z_IZ_O Z_OW_O.Female 0.25
ZW Z_IW_O Z_OZ_O Z_IZ_O.Male 0.50
ZW Z_IW_O Z_OZ_O Z_OW_O.Female 0.50
ZW Z_OW_I Z_IZ_I Z_IZ_O.Male 0.50
ZW Z_OW_I Z_IZ_I Z_IW_I.Female 0.50
ZW Z_OW_I Z_IZ_O Z_IZ_O.Male 0.25
ZW Z_OW_I Z_IZ_O Z_OZ_O.Male 0.25
ZW Z_OW_I Z_IZ_O Z_IW_I.Female 0.25
ZW Z_OW_I Z_IZ_O Z_OW_I.Female 0.25
ZW Z_OW_I Z_OZ_O Z_OZ_O.Male 0.50
ZW Z_OW_I Z_OZ_O Z_OW_I.Female 0.50
ZW Z_OW_O Z_IZ_I Z_IZ_O.Male 0.50
ZW Z_OW_O Z_IZ_I Z_IW_O.Female 0.50
ZW Z_OW_O Z_IZ_O Z_IZ_O.Male 0.25
ZW Z_OW_O Z_IZ_O Z_OZ_O.Male 0.25
ZW Z_OW_O Z_IZ_O Z_IW_O.Female 0.25
ZW Z_OW_O Z_IZ_O Z_OW_O.Female 0.25
ZW Z_OW_O Z_OZ_O Z_OZ_O.Male 0.50
ZW Z_OW_O Z_OZ_O Z_OW_O.Female 0.50
Code
make_mating_table(gene_location = "cytoplasmic") %>% 
  select(locus_type, everything()) %>% 
  rename(zygote_type = type) %>% 
  kable() %>% 
  kable_styling() %>% 
  scroll_box(height = "500px")
locus_type Female_genotype Male_genotype zygote_type zygote_freq
cytoplasmic C_I C_I C_I.Female 0.5
cytoplasmic C_I C_I C_I.Male 0.5
cytoplasmic C_I C_O C_I.Female 0.5
cytoplasmic C_I C_O C_I.Male 0.5
cytoplasmic C_O C_I C_O.Female 0.5
cytoplasmic C_O C_I C_O.Male 0.5
cytoplasmic C_O C_O C_O.Female 0.5
cytoplasmic C_O C_O C_O.Male 0.5

Load the inheritance schemes, to speed up the simulation

Code
# load the possible inheritance schemes

offspring_genotypes_autosome <- 
  make_mating_table(gene_location = "autosomal") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  as.data.table()

offspring_genotypes_X <- 
  make_mating_table(gene_location = "XY") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  filter(!str_detect(Male_genotype, "Y_I")) %>% 
  as.data.table()

offspring_genotypes_Y <- 
  make_mating_table(gene_location = "XY") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  filter(!str_detect(Female_genotype, "X_I"),
         !str_detect(Male_genotype, "X_I")) %>% 
  as.data.table()

offspring_genotypes_Z <- 
  make_mating_table(gene_location = "ZW") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  filter(!str_detect(Female_genotype, "W_I")) %>% 
  as.data.table()

offspring_genotypes_W <- 
  make_mating_table(gene_location = "ZW") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  filter(!str_detect(Female_genotype, "Z_I"),
         !str_detect(Male_genotype, "Z_I")) %>% 
  as.data.table()

offspring_genotypes_cytoplasmic <- 
  make_mating_table(gene_location = "cytoplasmic") %>% 
  select(1:4) %>% 
  rename(zygote_type = type) %>% 
  separate_wider_delim(zygote_type, names = c("Genotype", "Sex"), delim = ".") %>% 
  mutate(Sex = if_else(Sex == "Female", 1, 0)) %>% 
  as.data.table()

Build the main simulation function

Code
continuous_time_simulation <- function(row,
                                       parameters,
                                       inheritance_scheme){
  # package loading need to be done within the function to run on a cluster - local environment is not inherited by cluster
  library(data.table) 
  library(tidyverse)
  
  sample_vec <- function(x, ...) x[sample(length(x), ...)] # so we can sample from vectors with length 1 without this being interpreted as an integer
  
  print(paste("Doing row", row)) # this shows which row in the parameter space is being modelled
  
   prop_i_table <- data.table(time = numeric(), 
                             proportion_I = numeric(),
                            population_size = numeric()) # filled in as sim progresses
  
  keep_going <- TRUE # if the inbreeding allele fixes or goes extinct, this will change to false and the while loop will quit early
  
  Starting_pop_size <- parameters$Starting_pop_size[row] 
  N <- parameters$N[row] # constant
  number_mutants <- parameters$number_mutants[row] # constant at 1 
  baseline_mean_lifespan <- parameters$baseline_mean_lifespan[row] # constant at 2
  time_end <- parameters$time_end[row] # a cut-off point for each run 
  sex_expressed <- parameters$sex_expressed[row]
  chromosome <- parameters$chromosome[row]
  heterozygous_genotype <- parameters$heterozygous_genotype[row]
  homozygous_genotype <- parameters$homozygous_genotype[row]
  hemizygous_genotype <- parameters$hemizygous_genotype[row]
  #C <- parameters$C[row]
  v <- parameters$v[row]
  refractory_period <- parameters$refractory_period[row]
  D <- parameters$D[row]
  dominance <- parameters$dominance[row]
  parameter_space_ID <- parameters$parameter_space_ID[row]
  
  # define the starting genotypes for each sex so the population table can be built
  
  Female_starting_genotype <- inheritance_scheme[.N]$Female_genotype
  
  Male_starting_genotype <- inheritance_scheme[.N]$Male_genotype
  
  # Set the number of breeding sites
  
  breeding_sites <- 0.5*Starting_pop_size
  
  # Initialize the Individual_ID  and Family_ID counters
  
  Individual_ID_counter <- Starting_pop_size
  
  Family_ID_counter <- Starting_pop_size/N # family size equals the no. offspring produced by a single female
  
  # the simulation tracks the population via a data.table
  
  # create the starting population - note that females are sex = 1 and males are sex = 0
  
  population <-
    data.table(mortality_rate = 1/baseline_mean_lifespan,
               Sex = rbinom(n = Starting_pop_size, 1, prob = 0.5),
               birth_time = 0,
               matings = 0,
               reproduced = 0,
               mated_with = "NA",
               inbred_mating = 0,
               refractory_period_end = 0,
               relative_encountered = 0
    )[, `:=` (Genotype = ifelse(Sex > 0, Female_starting_genotype, Male_starting_genotype),
              breeding = Sex,
              Individual_ID = .I,
              Family_ID = rep(1:(.N/N), each = N, length.out = .N))]
  
  # seed population with the inbreeding allele
  
  # original - random alleles are mutated 
  
  #population[sample(which(str_detect(Genotype, pattern = chromosome)), size = number_mutants, replace = F),
   #         Genotype := str_replace(Genotype, pattern = paste0(chromosome, "_O"), replacement = paste0(chromosome, "_I"))]
  
  # current - entire family is mutated on one chromosome, unless the inbreeding allele is on a hemizygous chromosome, then the minimum number of families are mutated to make 5 mutations
  
  if(chromosome != "W" & chromosome != "Y"){
    population[Family_ID %in% sample(unique(population[, Family_ID]), size = 1),
               Genotype := str_replace(Genotype, pattern = paste0(chromosome, "_O"), replacement = paste0(chromosome, "_I"))]
  } else{
    # number_mutants is loaded from the parameter table
    mutations_made <- 0
    
    # Iterate through unique families as needed
    Families <- sample(unique(population$Family_ID), Starting_pop_size/N)
    Family_index <- 1
    
    while (mutations_made < number_mutants && Family_index <= length(Families)) {
      current_family <- Families[Family_index]
      
      # Identify mutate-able rows from the current family
      Family_rows <- population[Family_ID == current_family & str_detect(Genotype, pattern = chromosome)]
      
      # Calculate how many mutations we can make in this family
      changes_needed <- number_mutants - mutations_made
      changes_to_make <- min(nrow(Family_rows), changes_needed)
      
      # Select the individuals to mutate
      rows_to_modify <- Family_rows[1:changes_to_make]
      
      # Modify the genotype column in the selected rows
      population[rows_to_modify, 
                 Genotype := str_replace(Genotype, pattern = paste0(chromosome, "_O"), replacement = paste0(chromosome, "_I")),
                 on = "Individual_ID"]
      
      # Update the count of changes made
      mutations_made <- nrow(population[str_detect(Genotype, "I")])
      
      # Move to the next group
      Family_index <- Family_index + 1
    } 
  }
  
  # Determining the next event
  
  # check when the next death occurs: this is the sum of the mortality rates for all individuals in the population
  
  next_death <- rexp(n = 1, rate = sum(population[, mortality_rate]))
  
  # check when the next receptive female-male encounter occurs
  
  receptive_females <- population[Sex > 0]$Individual_ID # everyone's receptive at the start
  receptive_males <- population[Sex < 1]$Individual_ID
  
  # below is the slow way, which I still include to show it's equivalent to the fast way
  #encounter_possibilities <- 
  # CJ(Female_ID = population[Sex > 0]$Individual_ID,
  #   Male_ID = population[Sex < 1]$Individual_ID)[, encounter_rate := v/number_females]
  
  # Find the time the next encounter occurs: plug the sum of the rates into the exponential function. 
  # The population level encounter rate is the product of the rate at which a single male finds a single female, the number of receptive females in the population, and the number of receptive males in the population

  # 
  
  next_encounter <- rexp(n = 1, rate = length(receptive_females)*length(receptive_males)*v)
  
  # Initialize the timer t to the first encounter
  
  t <- pmin(next_death, next_encounter)
  
  # With the initial population ready to go and the first event found, start the timer and let the simulation run. In short, time progresses as events occur. Events can trigger state changes for the individuals in the population, leading to death, mating and offspring production.
  
  while (t <= time_end & keep_going){
    
    # death and mating: what type of encounter happens at time t
    
    if(next_death < next_encounter){
      
      who_died <- population[sample(.N, 1, prob = mortality_rate)]$Individual_ID 
      
      # remove individual from population table
      
      population <- population[Individual_ID != who_died]
      
      # check if mortality event frees up breeding site
      
      current_breeders <- sum(population$breeding > 0)
      
      # If there is an available breeding site, and at least one female to fill it, recruit a new breeder
      
      if(current_breeders < breeding_sites && sum(population$Sex > 0 & population$breeding < 1) > 0){
        
        # assign the new breeders
        
        population <- population[sample_vec(which(breeding < 1 & Sex > 0),
                                        size = 1), # note that all living females have equal prob of becoming a breeder
                                 breeding := 1]
      }
    } 
    else{
      
      which_encounter <- c(sample_vec(receptive_females, 1), sample_vec(receptive_males, 1))
      
      mates <- NULL # reset this every time as a safeguard
      
      # Opposite sex encounters
      
      female <- population[Individual_ID %chin% which_encounter[1]]
      male <- population[Individual_ID %chin% which_encounter[2]]
      
      # First determine if a homogametic individual, heterozygous for the I allele, will inbreed on this occasion.  
      
      heterozygote_inbreeds <- rbinom(1, 1, prob = dominance)
      
      # find a brother for the female to encounter, if required
      
      brother_encountered <-
        population[Family_ID %chin% female$Family_ID
                   # assign a receptive sibling - this is the behavioural trait under selection: 
                   # should an individual accept an inbred mating early in life occur or not.
        ][Sex < 1 & t > refractory_period_end, .SD[sample(.N, 1)]] 
      
      if(nrow(brother_encountered)>0 &
         female$relative_encountered < 1){
        
        # do inbreeding
        
        if(# female heterozygotes
          sex_expressed > 0 & 
          str_detect(female$Genotype, heterozygous_genotype) & 
          heterozygote_inbreeds > 0 |
          # female homozygotes
          sex_expressed > 0 & 
          str_detect(female$Genotype, homozygous_genotype) |
          # female hemizygotes
          sex_expressed > 0 & 
          str_detect(female$Genotype, hemizygous_genotype) |
          # male heterozygotes
          sex_expressed < 1 & 
          str_detect(brother_encountered$Genotype, heterozygous_genotype) & 
          heterozygote_inbreeds > 0 |
          # male homozygotes
          sex_expressed < 1 & 
          str_detect(brother_encountered$Genotype, homozygous_genotype) |
          # male hemizygotes
          sex_expressed < 1 &
          str_detect(brother_encountered$Genotype, hemizygous_genotype)){
          
          # if mating occurred, update the population
          
          mates <- rbindlist(list(female, brother_encountered))
          
          population[mates,
                     `:=`(matings = matings + 1,
                          inbred_mating = fifelse(Sex > 0, 1, 0),
                          mated_with = fifelse(Sex > 0, brother_encountered$Genotype, NA),
                          refractory_period_end = fifelse(Sex < 1, t + (refractory_period * baseline_mean_lifespan), 0)),
                     on = .(Individual_ID)]
        }
        
        # females that had not encountered a relative early in life are coded to have now done so
        
        population[female,
                   relative_encountered := 1,
                   on = .(Individual_ID)]
        
      }
      
      if(female$relative_encountered < 1 &
         nrow(brother_encountered) < 1){
      
       # females that had no receptive brother to encounter are recorded as having had their chance for inbreeding early in life. When the male refractory period != 0, this is unlikely (because all siblings are produced at the same time) but possible. Most commonly, this will occur when a female produces an all-female brood (0.03125 probability)
        
        population[female,
                   relative_encountered := 1,
                   on = .(Individual_ID)]
      }
      
      # If the individual has already encountered their sibling, don't swap and let encounter proceed. 
      
      # If the pair happen to be from the same family, assume individuals can recognise this and depending on genotype, avoid or accept the inbred mating
      
      if(# female heterozygotes
        sex_expressed > 0 & # female expressed allele
        female$relative_encountered > 0 & # already encountered their relative
        female$Family_ID == male$Family_ID & # they've run into another relative 
        str_detect(female$Genotype, heterozygous_genotype) & # they carry 1 inbreeding allele
        heterozygote_inbreeds > 0 | # the allele is expressed
        # female homozygotes
        sex_expressed > 0 & # female expressed allele
        male$relative_encountered > 0 & # already encountered relative
        female$Family_ID == male$Family_ID & # # they've run into another relative
        str_detect(female$Genotype, homozygous_genotype) |
        # female hemizygotes
        sex_expressed > 0 & 
        male$relative_encountered > 0 & # already encountered relative
        female$Family_ID == male$Family_ID & # # they've run into another relative
        str_detect(female$Genotype, hemizygous_genotype) |
        # male heterozygotes
        sex_expressed < 1 & # male expressed allele
        female$relative_encountered > 0 & # already encountered relative
        female$Family_ID == male$Family_ID & # they've encountered another relative 
        str_detect(male$Genotype, heterozygous_genotype) & # they carry 1 inbreeding allele
        heterozygote_inbreeds > 0 | # the allele is expressed
        # male homozygotes
        sex_expressed < 1 & # male expressed allele
        female$relative_encountered > 0 & # already encountered their relative
        female$Family_ID == male$Family_ID & # # they've run into another relative
        str_detect(male$Genotype, homozygous_genotype) |
        # male hemizygotes
        sex_expressed < 1 & # male expressed allele
        female$relative_encountered > 0 & # already encountered their relative
        female$Family_ID == male$Family_ID & # # they've run into another relative
        str_detect(male$Genotype, hemizygous_genotype)){
        
        mates <- rbindlist(list(female, male))
        
        population[mates,
                   `:=`(matings = matings + 1,
                        inbred_mating = fifelse(Sex > 0, 1, 0),
                        mated_with = fifelse(Sex > 0, male$Genotype, NA),
                        refractory_period_end = fifelse(Sex < 1, t + (refractory_period * baseline_mean_lifespan), 0)),
                   on = .(Individual_ID)]
      }
      
      # encounters between non-relatives
      
      if(female$Family_ID != male$Family_ID &
         female$relative_encountered > 0){
        
        mates <- rbindlist(list(female, male))
        
        population[mates,
                   `:=`(matings = matings + 1,
                        mated_with = fifelse(Sex > 0, male$Genotype, NA),
                        refractory_period_end = fifelse(Sex < 1, t + (refractory_period * baseline_mean_lifespan), 0)),
                   on = .(Individual_ID)]
      }
      
    }
    
    # Are there consequences of death and mating: reproduction
    
    # check if a female can now produce offspring, either because they're previously mated and have secured a breeding site or because they already held a breeding site and have now mated
    
    new_mated_breeder <- population[Sex > 0 & matings > 0 & breeding > 0 & reproduced < 1, 
                                    .(Individual_ID,
                                      inbred_mating,
                                      Female_genotype = Genotype,
                                      Male_genotype = mated_with)]
    
    if(nrow(new_mated_breeder) > 0) {
      # add offspring to the population. Each mated female that holds a breeding site produces N offspring
      offspring <- 
        new_mated_breeder[inheritance_scheme, 
                          on = .(Female_genotype = Female_genotype,
                                 Male_genotype = Male_genotype), 
                          nomatch = NULL, allow.cartesian  = TRUE
        ][, .SD[sample(.N,
                       size = N, 
                       prob = zygote_freq, 
                       replace = T)]
        ][, Family_ID := .GRP + Family_ID_counter # assign these offspring to a new family 
        ][, .(Genotype, 
              Sex, 
              inbred_mating,
              Family_ID)
        ][, `:=`(mortality_rate = fifelse(inbred_mating > 0, 
                                         1/(baseline_mean_lifespan + D), # the cost of inbreeding: D <= 0
                                         1/baseline_mean_lifespan), # outbred offspring mortality rate
                 birth_time = t,
                 breeding = 0,
                 matings = 0,
                 reproduced = 0,
                 mated_with = "NA",
                 refractory_period_end = t,
                 relative_encountered = 0,
                 Individual_ID = .I + Individual_ID_counter,
                 inbred_mating = 0)]
      
      # bind the offspring table to the existing population table and update which females have reproduced 
      
      population <- rbindlist(list(population, offspring), use.names = TRUE
      )[new_mated_breeder, reproduced := 1, on = .(Individual_ID)]
      
      # update the Individual_ID counter
      Individual_ID_counter <- max(population$Individual_ID)
      Family_ID_counter <- max(population$Family_ID) 
    }
    
    # Calculate the frequency of the I allele, quit early if I fixes or goes extinct
    
    # calc allele freq if autosomal locus   
    if(chromosome == "A"){
      prop_i <-
        (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) + 
           2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]))/ (nrow(population)*2)
    }
    
    # calc allele freq if hemizygous locus: W, Y or cytoplasmic
    if(chromosome == "W" | chromosome == "Y" | chromosome == "C"){
      prop_i <-
        (length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)])/ 
           length(population$Genotype[str_detect(population$Genotype, chromosome)]))
    }
    
    # calc allele freq if diploid in one sex and haploid in the other: X and Z
    if(chromosome == "X" | chromosome == "Z"){
      prop_i <-
        if(hemizygous_genotype == "X_IY_O"){
          (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) + 
             2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]) +
             length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)]))/ 
            (nrow(population[Sex > 0])*2 + nrow(population[Sex < 1]))}
      else{
        (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) + 
           2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]) +
           length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)]))/ 
          (nrow(population[Sex < 1])*2 + nrow(population[Sex > 0]))}
    }
    
    # this is a diagnostic to make sure the model is running well - it can be commented out when running the big simulation
    prop_i_table <- rbindlist(list(prop_i_table, list(t, prop_i, nrow(population))))
    
    print(paste0("Population size = ", nrow(population), ", breeders = ", sum(population$breeding > 0), 
                 ", time = ", round(t, 3)))
    
    if(prop_i > 0.5 | prop_i < 0.0001 | nrow(population) < 2) keep_going <- FALSE
    
    # Move t to next encounter
    
    # determining the next event
    
    # check when the next death occurs
    
    next_death <- rexp(n = 1, rate = sum(population[, mortality_rate]))
    
    # check when the next receptive female-male encounter occurs
    
    receptive_females <- population[Sex > 0 & matings < 1]$Individual_ID
    receptive_males <- population[Sex < 1 & refractory_period_end <= t]$Individual_ID
    
    if(length(receptive_females)*length(receptive_males) > 0){
      
      # check the time: this is the sum of the rates, because each male finds each female at the same rate
      
      next_encounter <- rexp(n = 1, rate = length(receptive_females)*length(receptive_males)*v)
      
      # Initialize the timer t to the next encounter
      
    }else{next_encounter <- 10^3} # make this a number that will never be exceeded by a death time
    
    t <- t + pmin(next_death, next_encounter) 
    
  }
  finish_time <- t
  final_pop_size = nrow(population)

# cut the i table down to 1 row per 0.1 time increment, join with the relevant parameter space info and save as a csv.  
    
prop_i_table %>%
  mutate(time_group = floor(time / 0.1) * 0.1) %>%
  group_by(time_group) %>%
  slice(1) %>%
  ungroup() %>%
  select(-time_group) %>%
  bind_rows(prop_i_table %>% filter(row_number()==n())) %>%
  distinct() %>% 
  bind_cols(parameters[row, ] %>% 
    select(-c(heterozygous_genotype, homozygous_genotype, hemizygous_genotype,
              number_mutants, baseline_mean_lifespan, N, time_end))) %>% 
    write_csv(paste("sim_results/rowID_", parameter_space_ID, chromosome, ".csv", sep = ""))
}

Define the parameter space

The effect of the relatedness coefficient on the fitness of an inbreeding allele when there is non-zero inbreeding depression is well described. Therefore, we choose to set \(r = 0.5\) for all inbreeding events and instead focus on the effect of varying 1) \(v\), 2) \(\delta\), 3) the male cost to mating, here coded as a refractory period following mating, and 4) the chromosome upon which the the allele is found.

Populations also differ in the number of breeding females they can support, in order to standardise the number of chromosomes carrying the focal locus. For example, when the inbreeding locus is found on an autosome, a population of 200 individuals harbours 400 copies of the inbreeding locus. If the allele is instead found on a cytoplasmic chromosome (and we assume no heteroplasmy), each individual only carries one copy of the locus and a population size of 400 individuals is needed to equalise the population size of chromosomes. \(v\) is adjusted in accordance with the individual population size, so as to ensure that males meet the same number of individuals on average across simulation runs - this is akin to keeping the population density constant between simulation runs.

Code
resolution <- 30
starting_pop_size_autosomes <- 200 # both sexes harbour two copies of each autosomal chromosome = 400 autosomal haplotypes

parameters <-
  expand_grid(
    #C = c(2, 10),
    chromosome = c("A", "X", "Y", "Z", "W", "C"),
    v = c(0.02, 0.05, 2),
    D = seq(0, -0.99, length = resolution), # inbreeding depression
    refractory_period = seq(0, 1, length = 15)
  ) %>% 
  full_join(tibble(chromosome = c("A", "A", "A", "A", "A", "A",
                                  "X", "X", "X", "X",
                                  "Y",
                                  "Z", "Z", "Z", "Z",
                                  "W",
                                  "C", "C"),
                   sex_expressed = c(0, 0, 0, 1, 1, 1,
                                     0, 1, 1, 1,
                                     0,
                                     0, 0, 0, 1,
                                     1,
                                     0, 1),
                   dominance = c(0, 0.5, 1, 0, 0.5, 1,
                                 1, 0, 0.5, 1, 
                                 1, 
                                 0, 0.5, 1, 1,
                                 1,
                                 1, 1)) %>% 
              mutate(heterozygous_genotype = case_when(chromosome == "A" ~ "A_IA_O",
                                                       chromosome == "X" ~ "X_IX_O",
                                                       chromosome == "Y" ~ "NA",
                                                       chromosome == "Z" ~ "Z_IZ_O",
                                                       chromosome == "W" ~ "NA",
                                                       chromosome == "C" ~ "NA"),
                     homozygous_genotype = case_when(chromosome == "A" ~ "A_IA_I",
                                                     chromosome == "X" ~ "X_IX_I",
                                                     chromosome == "Y" ~ "NA",
                                                     chromosome == "Z" ~ "Z_IZ_I",
                                                     chromosome == "W" ~ "NA",
                                                     chromosome == "C" ~ "NA"),
                     hemizygous_genotype = case_when(chromosome == "A" ~ "NA",
                                                     chromosome == "X" ~ "X_IY_O",
                                                     chromosome == "Y" ~ "X_OY_I",
                                                     chromosome == "Z" ~ "Z_IW_O",
                                                     chromosome == "W" ~ "Z_OW_I",
                                                     chromosome == "C" ~ "C_I"),
                     Starting_pop_size = case_when(chromosome == "A" ~ starting_pop_size_autosomes,
                                                   chromosome == "X" | chromosome == "Z" ~ 
                                                     starting_pop_size_autosomes + 0.33*starting_pop_size_autosomes,
                                                   chromosome == "Y" | chromosome == "W" ~ starting_pop_size_autosomes*4,
                                                   chromosome == "C" ~ starting_pop_size_autosomes*2)),
            relationship = "many-to-many", by = "chromosome") %>% 
  mutate(baseline_mean_lifespan = 1,
         v = v * (starting_pop_size_autosomes / Starting_pop_size),
         N = 5, # subject to change
         number_mutants = N, # this is the size of a family
         time_end = 1000, # with avg lifespan = 1, this is ~ roughly 1000 gens
         parameter_space_ID = row_number())

parameters_autosome <- parameters %>% filter(chromosome == "A")
parameters_X <- parameters %>% filter(chromosome == "X")
parameters_Y <- parameters %>% filter(chromosome == "Y")
parameters_Z <- parameters %>% filter(chromosome == "Z")
parameters_W <- parameters %>% filter(chromosome == "W")
parameters_C <- parameters %>% filter(chromosome == "C")

Run the simulation

We run one simulation per row in the parameter space table.

Code
parameters_autosome <- parameters_autosome %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_X <- parameters %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_Z <- parameters_Z %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_Y <- parameters_Y %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_W <- parameters_W %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs
parameters_C <- parameters_C %>% slice_sample(prop = 1) # shuffle to equalise workload across jobs

# run the simulation in parallel using all local CPUs on your computer. The code below uses the autosomal parameter space as an example.

if(!file.exists("results/complete_results.csv")){
  
  # get the number of cores on the computer 
  n.cores <- parallel::detectCores()
  
  # make the cluster
  cluster <- parallel::makeCluster(n.cores)
  
  # split parameters across cluster
  parallel::clusterExport(cluster, paste("parameters_Y"))
  
  parallel::parLapply(cluster,
                      1:nrow(parameters_Y),
                      continuous_time_simulation,
                      parameters_Y,
                      offspring_genotypes_Y)
}

Load the results

Code
# build a function to load the individual runs and join them into a single tibble

load_results <- function(chromosome){
  
  files <-
    list.files(path = "sim_results") %>% 
    str_subset(chromosome)
  
  paste("sim_results/", files, sep = "") %>% 
    vroom() 
}

if(!file.exists("results/complete_results.csv")){
  cytoplasmic_results <- load_results("C")
  y_sim_results <- load_results("Y") 
  X_sim_results <- load_results("X")
  Z_sim_results <- load_results("Z")
  W_sim_results <- load_results("W")
  A_sim_results <- load_results("A")
  
  sim_results <-  
    bind_rows(A_sim_results,
              y_sim_results,
              X_sim_results,
              Z_sim_results,
              W_sim_results,
              cytoplasmic_results) 
  
  sim_results %>% 
    write_csv("results/complete_results.csv")} else{
      sim_results <- read_csv("results/complete_results.csv")
    }

Plot the results

Simulation diagnostics

Code
# load a palette used in all tabs
temp <- pnw_palette("Shuksan2",100)
Code
A_allele_tracking <-
  sim_results %>%
  filter(chromosome == "A") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
                                   v = round(v, 2)) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200), breaks = c(0, 50, 100, 150)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

A_pop_tracking <-
  sim_results %>%
  filter(chromosome == "A") %>%   
    mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
                                   v = round(v, 2)) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
A_pop_tracking / A_allele_tracking

Code
X_allele_tracking <-
  sim_results %>%
  filter(chromosome == "X") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200), breaks = c(0, 50, 100, 150)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(~sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

X_pop_tracking <-
  sim_results %>%
  filter(chromosome == "X") %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
X_pop_tracking/X_allele_tracking

Code
Z_allele_tracking <-
  sim_results %>%
  filter(chromosome == "Z") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
                                   v = round(v, 2)) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200), breaks = c(0, 50, 100, 150)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

Z_pop_tracking <-
  sim_results %>%
  filter(chromosome == "Z") %>%   
    mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
                                   v = round(v, 2)) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
Z_pop_tracking /Z_allele_tracking

Code
Y_allele_tracking <-
  sim_results %>%
  filter(chromosome == "Y") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 200), breaks = c(0, 50, 100, 150)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

Y_pop_tracking <-
  sim_results %>%
  filter(chromosome == "Y") %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
Y_pop_tracking/Y_allele_tracking

Code
W_allele_tracking <-
  sim_results %>%
  filter(chromosome == "W") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 125), breaks = c(0, 25, 50, 75, 100)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

W_pop_tracking <-
  sim_results %>%
  filter(chromosome == "W") %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
W_pop_tracking/W_allele_tracking

Code
C_allele_tracking <-
  sim_results %>%
  filter(chromosome == "C") %>% 
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = proportion_I)) + 
  geom_line(aes(group = parameter_space_ID, colour = D_prop), alpha = 0.4) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 125), breaks = c(0, 25, 50, 75, 100)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Frequency of _I_ allele",
       colour = "Inbreeding depression") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))

C_pop_tracking <-
  sim_results %>%
  filter(chromosome == "C") %>%   
  mutate(D_prop = D*-1,
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females")) %>% 
  ggplot(aes(x = time, y = population_size)) + 
  geom_line(aes(group = parameter_space_ID, colour = refractory_period), alpha = 0.6) +
  #geom_smooth(linewidth = 2) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) + 
  #scale_y_continuous(expand = c(0, 0), limits = c(0, 0.5)) +
  scale_colour_gradientn(colours = temp) +
  labs(x = "Time",
       y = "Population size",
       colour = "Mating system") +
  facet_wrap(sex_expressed~v, nrow = 2,
             labeller = 
               label_glue('Sex expressing inbreeding allele: {`sex_expressed`}\nMale search efficiency: {`v`}')) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.title.y = element_markdown(),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        axis.text = element_text(size = 12))
        
C_pop_tracking/C_allele_tracking

When is inbreeding preference selected?

Code
data <-
  sim_results %>%
  filter(chromosome == "A") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.02 ~ "Very low",
                           v == 0.05 ~ "Low",
                           .default = "High"))

A_heatmap <-
  data %>%
  #filter(dominance == 1) %>% 
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 0.7) + 
  geom_vline(data = data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5*1)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*1)/(1 + 0.5*1)), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], 
                               pnw_palette("Shuksan2", n = 5)[4], 
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "Autosomal inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

A_heatmap

Code
data <-
  sim_results %>%
  filter(chromosome == "X") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         v = round(v, 3),
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.015 ~ "Very low",
                           v == 0.038 ~ "Low",
                           .default = "High"))

X_heatmap <-
  data %>%
  #filter(dominance == 1) %>% 
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 0.7) + 
  geom_vline(data = data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + 0.5*1)), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*0.5)/(1 + 0.5*0.5)), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], 
                               pnw_palette("Shuksan2", n = 5)[4], 
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "X-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

X_heatmap

Code
data <-
  sim_results %>%
  filter(chromosome == "Z") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         v = round(v, 3),
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.015 ~ "Very low",
                           v == 0.038 ~ "Low",
                           .default = "High"))

Z_heatmap <-
  data %>%
  #filter(dominance == 1) %>% 
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 0.7) + 
  geom_vline(data = data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1/(1 + (0.5*0.5))), linetype = 2, colour = "white", linewidth = .8) +
  geom_vline(data = data %>% filter(sex_expressed == "Females"),
             aes(xintercept = (0.5*1)/(1 + (0.5*1))), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], 
                               pnw_palette("Shuksan2", n = 5)[4], 
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "Z-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Z_heatmap

Code
data <-
  sim_results %>%
  filter(chromosome == "Y") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.005 ~ "Very low",
                           v == 0.0125 ~ "Low",
                           .default = "High"))

Y_heatmap <-
  data %>%
  #filter(dominance == 1) %>% 
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(data = data %>% filter(sex_expressed == "Males"),
             aes(xintercept = 1), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], 
                               pnw_palette("Shuksan2", n = 5)[4], 
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "Y-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

Y_heatmap

Code
data <-
  sim_results %>%
  filter(chromosome == "W") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>% 
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.005 ~ "Very low",
                           v == 0.0125 ~ "Low",
                           .default = "High"))

W_heatmap <-
  data %>%
  #filter(dominance == 1) %>% 
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(aes(xintercept = 0), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], 
                               pnw_palette("Shuksan2", n = 5)[4], 
                               pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "W-linked inbreeding alleles") +
  facet_wrap(sex_expressed~ v_cat, nrow = 1,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

W_heatmap

Code
data <- 
  sim_results %>%
  filter(chromosome == "C") %>%  
  group_by(parameter_space_ID) %>% 
  slice_tail() %>%
  ungroup() %>% 
  mutate(D_prop = D*-1,
         Fate = case_when(population_size < Starting_pop_size*0.1 ~ "Extinction",
                          proportion_I  > 0.5 ~ "Invades",
                          proportion_I  < 0.001 ~ "Purged"),
         sex_expressed = case_when(sex_expressed == 0 ~ "Males",
                                   .default = "Females"),
         v_cat = case_when(v == 0.01 ~ "Very low",
                           v == 0.025 ~ "Low",
                           .default = "High"))

C_heatmap <-
  data %>%
  ggplot(aes(x = D_prop, y = refractory_period)) +
  geom_blank() +
  geom_tile(aes(fill = Fate), alpha = 1) + 
  geom_vline(aes(xintercept = 0), linetype = 2, colour = "white", linewidth = .8) +
  scale_fill_manual(values = c(pnw_palette("Shuksan2", n = 5)[2], pnw_palette("Shuksan2", n = 5)[4], pnw_palette("Shuksan2", n = 5)[5])) +
  labs(x = ~delta~'(inbreeding depression)',
       y = "Male refractory period",
       fill = "Fate of _I_ allele",
       title = "Cytoplasmic inbreeding alleles") +
  facet_wrap(sex_expressed ~ v_cat, nrow = 2,
             labeller = label_glue('Sex: {`sex_expressed`}\nMale search efficiency: {`v_cat`}')) +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_y_continuous(expand = c(0, 0)) + 
  theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
        strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
        strip.text = element_text(size = 10),
        plot.title = element_text(hjust = 0.5, size = 16),
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        legend.title = element_markdown())

C_heatmap

Session information

Code
sessionInfo() %>% pander

R version 4.4.0 (2024-04-24)

Platform: x86_64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bench(v.1.1.3), profvis(v.0.3.8), ggnewscale(v.0.4.10), patchwork(v.1.2.0), ggtext(v.0.1.2), geomtextpath(v.0.1.4), stickylabeller(v.0.0.0.9100), tidybayes(v.3.0.6), PNWColors(v.0.1.0), wesanderson(v.0.3.7), MoMAColors(v.0.0.0.9000), MetBrewer(v.0.2.0), data.table(v.1.15.4), kableExtra(v.1.4.0), pander(v.0.6.5), vroom(v.1.6.5), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): svUnit(v.1.0.6), tidyselect(v.1.2.1), viridisLite(v.0.4.2), farver(v.2.1.2), fastmap(v.1.2.0), tensorA(v.0.36.2.1), digest(v.0.6.35), timechange(v.0.3.0), lifecycle(v.1.0.4), magrittr(v.2.0.3), posterior(v.1.5.0), compiler(v.4.4.0), rlang(v.1.1.4), tools(v.4.4.0), utf8(v.1.2.4), yaml(v.2.3.8), knitr(v.1.47), labeling(v.0.4.3), htmlwidgets(v.1.6.4), bit(v.4.0.5), xml2(v.1.3.6), abind(v.1.4-5), withr(v.3.0.0), grid(v.4.4.0), fansi(v.1.0.6), ggstream(v.0.1.0), colorspace(v.2.1-0), scales(v.1.3.0), isoband(v.0.2.7), cli(v.3.6.2), rmarkdown(v.2.27), crayon(v.1.5.2), generics(v.0.1.3), rstudioapi(v.0.16.0), tzdb(v.0.4.0), commonmark(v.1.9.1), parallel(v.4.4.0), vctrs(v.0.6.5), jsonlite(v.1.8.8), hms(v.1.1.3), arrayhelpers(v.1.1-0), bit64(v.4.0.5), systemfonts(v.1.1.0), ggdist(v.3.3.2), glue(v.1.7.0), cowplot(v.1.1.3), distributional(v.0.4.0), stringi(v.1.8.4), gtable(v.0.3.5), munsell(v.0.5.1), pillar(v.1.9.0), htmltools(v.0.5.8.1), R6(v.2.5.1), textshaping(v.0.4.0), evaluate(v.0.24.0), lattice(v.0.22-6), highr(v.0.11), markdown(v.1.13), backports(v.1.5.0), gridtext(v.0.1.5), Rcpp(v.1.0.12), svglite(v.2.1.3), coda(v.0.19-4.1), checkmate(v.2.3.1), xfun(v.0.44) and pkgconfig(v.2.0.3)